import sys
# Install gurobi
!{sys.executable} -m pip install -i https://pypi.gurobi.com gurobipy
!{sys.executable} -m pip install plotly geopy

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pandas as pd
from itertools import product

import gurobipy as gp
from gurobipy import GRB, quicksum

# These are ours
import problem
import visualization

Looking in indexes: https://pypi.gurobi.com
Requirement already satisfied: gurobipy in /opt/anaconda3/lib/python3.8/site-packages (9.1.1)
Requirement already satisfied: plotly in /opt/anaconda3/lib/python3.8/site-packages (4.14.3)
Requirement already satisfied: geopy in /opt/anaconda3/lib/python3.8/site-packages (2.1.0)
Requirement already satisfied: six in /opt/anaconda3/lib/python3.8/site-packages (from plotly) (1.15.0)
Requirement already satisfied: retrying>=1.3.3 in /opt/anaconda3/lib/python3.8/site-packages (from plotly) (1.3.3)
Requirement already satisfied: geographiclib<2,>=1.49 in /opt/anaconda3/lib/python3.8/site-packages (from geopy) (1.50)

def normalize(i, j):
    assert i != j
    return (i, j) if i < j else (j, i)

Problem configuration

from model import Autonomax, Config

# Run on the first |C| cities (mostly for testing)
C = 41

# What scenario we will be utilizing
scenario = 1
# The number of cities in the core net
NC = 3
# 1 if the core net should be a cycle; 0 if it shoul be a path.
Z = 1

autonomax = Autonomax(
    Config(
        cities=problem.cities[:C],
        distances=problem.D[:C, :C],
        demand=problem.B[scenario][:C],
        core_city_count=NC,
        core_net_is_cycle=Z,
    )
)

Academic license - for non-commercial use only - expires 2021-05-15
Using license file /Users/sjurwold/gurobi.lic

A = autonomax.model.getA()

count = lambda d: len(d) if isinstance(d, gp.tupledict) else 1

V = sum(count(v) for v in autonomax.variables.values())
C = sum(count(c) for c in autonomax.constraints.values())

print(f'V = {V}, C = {C}')

fig, ax = plt.subplots(1, figsize=(128, 128 * C/V))
ax.set_xticklabels([])
ax.set_yticklabels([])
plt.spy(A)


total = 0
print(f'\nCONSTRAINT SETS:')
for (name, constraint) in autonomax.constraints.items():
    constraints_in_set = count(constraint)
    print(f'{constraints_in_set:>5} | {name}')
    #plt.gcf().text(0.07, 0.69 - 0.4 * (total + 0.5 * count(constraint)) / C, name, fontsize=14)
    total += constraints_in_set
    plt.plot([0, V], [total - 0.5, total - 0.5], color='grey')
    
total = 0
print(f'\nVARIABLE SETS:')
for (name, variables) in autonomax.variables.items():
    variables_in_set = count(variables)
    print(f'{variables_in_set:>5} | {name}')
    #plt.gcf().text(0.1 + 0.8 * (total + 0.5*count(variables)) / V, 0.8, name, fontsize=14)
    total += count(variables)
    plt.plot([total - 0.5, total - 0.5], [0, C], color='grey')
    
plt.savefig('constraint-matrix.png', dpi=200)

V = 7667, C = 5292

CONSTRAINT SETS:
    1 | one_control_center
   41 | control_city_directly_connected
 1681 | reduce_control_center_symmetry
   41 | core_city_ub
    1 | cycle_or_path
   41 | disallow_core_tree
    1 | exactly_nc_core_cities
   41 | control_center_is_connected
 1640 | is_connectable
   82 | is_connected_timestep
   41 | connected_graph
   41 | conservation_of_flow
  820 | force_edge_if_flow
  820 | edge_cost_lb

VARIABLE SETS:
   41 | is_control_center
  820 | is_core_edge
   41 | is_core_city
  123 | is_connected_step
 3362 | is_connectable_step
  820 | flow
  820 | abs_flow
  820 | is_sub_edge
  820 | edge_cost

autonomax.model.Params.timeLimit = 600.0
autonomax.model.optimize()

Changed value of parameter timeLimit to 600.0
   Prev: inf  Min: 0.0  Max: inf  Default: inf
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 5292 rows, 7667 columns and 25297 nonzeros
Model fingerprint: 0x103a39e2
Model has 820 general constraints
Variable types: 2460 continuous, 5207 integer (5207 binary)
Coefficient statistics:
  Matrix range     [1e+00, 2e+05]
  Objective range  [1e+00, 2e+04]
  Bounds range     [1e+00, 1e+02]
  RHS range        [8e-01, 4e+01]
Presolve removed 89 rows and 1768 columns
Presolve time: 0.08s
Presolved: 5203 rows, 5899 columns, 25440 nonzeros
Variable types: 2460 continuous, 3439 integer (3439 binary)
Found heuristic solution: objective 78529.573156

Root relaxation: objective 1.176011e+03, 5395 iterations, 0.12 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 1176.01073    0  144 78529.5732 1176.01073  98.5%     -    0s
     0     0 2149.64079    0  207 78529.5732 2149.64079  97.3%     -    0s
H    0     0                    70223.874648 2149.64079  96.9%     -    0s
H    0     0                    44023.687886 2149.64079  95.1%     -    0s
     0     0 2477.24371    0  245 44023.6879 2477.24371  94.4%     -    0s
     0     0 2478.53134    0  272 44023.6879 2478.53134  94.4%     -    0s
     0     0 2478.63320    0  273 44023.6879 2478.63320  94.4%     -    0s
     0     0 2607.46916    0  328 44023.6879 2607.46916  94.1%     -    0s
     0     0 2614.29349    0  319 44023.6879 2614.29349  94.1%     -    0s
     0     0 2614.34693    0  322 44023.6879 2614.34693  94.1%     -    0s
H    0     0                    19375.073949 2614.34693  86.5%     -    1s
     0     0 2639.27518    0  321 19375.0739 2639.27518  86.4%     -    1s
     0     0 2643.53420    0  311 19375.0739 2643.53420  86.4%     -    1s
     0     0 2645.72876    0  320 19375.0739 2645.72876  86.3%     -    1s
     0     0 2650.80305    0  321 19375.0739 2650.80305  86.3%     -    1s
     0     0 2650.82315    0  321 19375.0739 2650.82315  86.3%     -    1s
     0     0 2668.19595    0  332 19375.0739 2668.19595  86.2%     -    1s
     0     0 2673.81690    0  336 19375.0739 2673.81690  86.2%     -    1s
     0     0 2673.98505    0  336 19375.0739 2673.98505  86.2%     -    1s
     0     0 2675.50290    0  321 19375.0739 2675.50290  86.2%     -    1s
     0     0 2676.08915    0  323 19375.0739 2676.08915  86.2%     -    1s
     0     0 2680.71297    0  298 19375.0739 2680.71297  86.2%     -    1s
     0     0 2681.64640    0  326 19375.0739 2681.64640  86.2%     -    1s
     0     0 2681.66345    0  322 19375.0739 2681.66345  86.2%     -    1s
     0     0 2684.84889    0  295 19375.0739 2684.84889  86.1%     -    1s
     0     0 2684.84889    0  235 19375.0739 2684.84889  86.1%     -    2s
H    0     0                    16829.684245 2685.37753  84.0%     -    2s
     0     2 2685.37753    0  233 16829.6842 2685.37753  84.0%     -    2s
H   18    13                    15666.371601 2842.05742  81.9%   480    4s
H   75    61                    15657.732090 2868.02137  81.7%   196    4s
H  125    96                    15657.732081 2868.02137  81.7%   137    4s
   132   106 3104.64701   15  277 15657.7321 2868.02137  81.7%   133    5s
H  270   200                    15114.042650 2868.02137  81.0%   114    6s
H  502   364                    15095.720047 2868.02137  81.0%  98.9    6s
H  503   364                    15069.227778 2868.02137  81.0%  99.1    6s
H  505   364                    15041.752975 2868.02137  80.9%  98.8    6s
H  661   481                    14069.040385 2868.02137  79.6%  93.9    7s
H  729   470                    13732.088544 2868.02137  79.1%  91.1    7s
  1452   987 9632.75470   47  235 13732.0885 2868.02137  79.1%  89.1   11s
  1475  1002 8554.82656   50  343 13732.0885 2868.02137  79.1%  87.7   15s
  1495  1016 8172.73998   28  395 13732.0885 2996.79994  78.2%  86.5   20s
  1516  1030 12081.3489   52  435 13732.0885 3083.76581  77.5%  85.3   25s
  1526  1037 3174.46733   12  408 13732.0885 3174.46733  76.9%   120   35s
  1533  1038 3230.06510   15  365 13732.0885 3199.20521  76.7%   136   41s
  1542  1043 3360.72560   17  334 13732.0885 3277.06727  76.1%   142   45s
  1646  1074 6366.56877   22  198 13732.0885 3512.53706  74.4%   163   50s
  1841  1138 7717.22693   27  193 13732.0885 3512.53706  74.4%   172   55s
  2166  1125 7130.29866   33  172 13732.0885 3512.53706  74.4%   185   60s
  2771  1068 4223.63351   28  210 13732.0885 3675.24044  73.2%   185   65s
  3026  1141 11357.3371   45   89 13732.0885 4015.42961  70.8%   189   70s
  3648  1286 8326.28011   34  158 13732.0885 4595.65379  66.5%   193   75s
H 3963  1312                    13726.786579 5080.74060  63.0%   192   77s
H 4122  1278                    13698.523104 5223.68602  61.9%   191   77s
H 4184  1212                    13698.523059 5223.68602  61.9%   191   77s
H 4225  1211                    13695.037175 5223.68602  61.9%   191   77s
H 4410  1345                    13695.037165 5326.77399  61.1%   192   79s
  4532  1423 7865.62664   42  104 13695.0372 5389.17046  60.6%   192   80s
H 4924  1498                    13695.037146 5815.22786  57.5%   192   82s
H 5020  1634                    13695.037129 5830.98039  57.4%   193   84s
  5257  1712 infeasible   42      13695.0371 5950.13765  56.6%   192   85s
H 5263  1712                    13695.037123 5950.13765  56.6%   192   85s
H 5329  1712                    13695.037118 6054.14969  55.8%   193   85s
H 5785  1918                    13695.037105 6329.09380  53.8%   193   87s
H 6068  1991                    13695.037095 6465.63348  52.8%   191   89s
H 6133  1991                    13695.037086 6541.83051  52.2%   191   89s
  6217  2013 7943.26840   34  127 13695.0371 6586.01184  51.9%   191   90s
H 6630  2125                    13695.037079 6865.08626  49.9%   195   93s
H 6646  2125                    13695.037074 6873.36168  49.8%   195   93s
  6852  2143 9630.26145   34  151 13695.0371 7169.02866  47.7%   199   96s
H 6868  2143                    13695.037071 7169.02866  47.7%   199   96s
  7470  2192 9056.23701   40  113 13695.0371 7601.95145  44.5%   204  100s
  7853  2132     cutoff   32      13695.0371 7888.07820  42.4%   207  128s
  8110  2046     cutoff   32      13695.0371 8098.29187  40.9%   210  185s
H 8688  1876                    13695.037067 8638.78282  36.9%   216  188s
H 8747  1876                    13695.037015 8673.61925  36.7%   216  188s
  8949  1711 10844.0192   42   96 13695.0370 8884.23798  35.1%   218  190s
  9322  1518     cutoff   53      13695.0370 9222.04478  32.7%   219  218s
H 9339  1518                    13695.037009 9222.04478  32.7%   219  218s
H 9678  1518                    13695.037003 9534.92139  30.4%   220  218s
  9722  1178     cutoff   47      13695.0370 9654.61621  29.5%   220  220s
H 9731  1178                    13695.036977 9654.61621  29.5%   220  220s
H 9797  1178                    13695.036928 9657.06722  29.5%   220  220s
H10227   741                    13695.036915 10209.8330  25.4%   219  222s
H10377   741                    13695.036867 10229.2922  25.3%   219  222s
H11094   249                    13695.036865 11216.1210  18.1%   213  224s
 11285   278 infeasible   36      13695.0369 12293.8951  10.2%   211  225s

Cutting planes:
  Gomory: 31
  Cover: 290
  MIR: 308
  Flow cover: 567
  Inf proof: 4
  Zero half: 3
  RLT: 115
  Relax-and-lift: 73

Explored 14892 nodes (2419021 simplex iterations) in 226.54 seconds
Thread count was 8 (of 8 available processors)

Solution count 10: 13695 13695 13695 ... 13695

Optimal solution found (tolerance 1.00e-04)
Best objective 1.369503686545e+04, best bound 1.369503686545e+04, gap 0.0000%

Solution

city_info = pd.DataFrame(autonomax.city_info())
city_info

	Index	Name	IsCoreCity	IsControlCenter	Demand	IngoingFlow	OutgoingFlow
0	0	Boden	False	False	4.4	2.700000e+00	7.100000
1	1	Borås	False	False	3.4	2.100000e+01	24.400000
2	2	Eskilstuna	False	False	3.2	6.200000e+00	9.400000
3	3	Falun	False	False	4.6	1.136868e-13	4.600000
4	4	Gävle	False	False	2.6	2.380000e+01	26.400000
5	5	Göteborg	False	False	6.4	1.207923e-13	6.400000
6	6	Halmstad	False	False	4.0	7.400000e+00	11.399999
7	7	Haparanda	False	False	1.1	0.000000e+00	1.100000
8	8	Helsingborg	False	False	3.0	4.400000e+00	7.400000
9	9	Hudiksvall	False	False	1.5	1.880000e+01	20.299999
10	10	Jönköping	False	False	2.2	2.440000e+01	26.600000
11	11	Kalmar	False	False	1.9	5.755396e-13	1.900000
12	12	Karlskrona	False	False	3.2	2.984279e-13	3.200000
13	13	Karlstad	False	False	4.4	8.100187e-13	4.400000
14	14	Kiruna	False	False	2.7	2.700062e-13	2.700000
15	15	Kristianstad	False	False	3.9	4.902745e-13	3.900000
16	16	Lidköping	False	False	3.4	2.700000e+00	6.100000
17	17	Linköping	True	False	4.6	1.161000e+02	120.700000
18	18	Luleå	False	False	0.9	8.200000e+00	9.100000
19	19	Malmö	False	False	0.8	3.600000e+00	4.400000
20	20	Motala	True	False	3.6	1.073000e+02	110.900000
21	21	Norrköping	True	True	3.9	1.704000e+02	53.600000
22	22	Nyköping	False	False	1.7	3.860000e+01	40.300000
23	23	Sandviken	False	False	3.5	9.592327e-13	3.500000
24	24	Skellefteå	False	False	1.8	9.100000e+00	10.899999
25	25	Skövde	False	False	1.0	6.100000e+00	7.100000
26	26	Stockholm	False	False	9.2	2.940000e+01	38.600000
27	27	Sundsvall	False	False	1.0	1.780000e+01	18.799999
28	28	Trelleborg	False	False	3.6	1.200817e-12	3.600000
29	29	Uddevalla	False	False	1.2	1.008971e-12	1.200000
30	30	Umeå	False	False	2.4	1.090000e+01	13.299999
31	31	Uppsala	False	False	3.0	2.640000e+01	29.400000
32	32	Varberg	False	False	3.2	1.140000e+01	14.600000
33	33	Vetlanda	False	False	2.8	1.060000e+01	13.400000
34	34	Vänersborg	False	False	1.5	1.200000e+00	2.700000
35	35	Västervik	False	False	3.3	1.900000e+00	5.200000
36	36	Västerås	False	False	1.6	4.600000e+00	6.200000
37	37	Växjö	False	False	3.5	7.100000e+00	10.600000
38	38	Örebro	False	False	2.2	4.400000e+00	6.600000
39	39	Örnsköldsvik	False	False	2.2	1.330000e+01	15.499999
40	40	Östersund	False	False	2.3	9.166001e-13	2.300000

edge_info = pd.DataFrame(autonomax.edge_info())
edge_info

	From	To	Type	Flow	Cost	Distance
0	Skellefteå	Umeå	SUB	10.899999	413.098267	111
1	Eskilstuna	Västerås	SUB	-6.200000	65.283353	43
2	Gävle	Hudiksvall	SUB	-20.299999	832.846872	118
3	Norrköping	Nyköping	SUB	-40.300000	572.920459	58
4	Nyköping	Stockholm	SUB	-38.600000	1052.199991	90
5	Linköping	Norrköping	CORE	120.700000	440.000000	44
6	Karlskrona	Växjö	SUB	3.200000	114.243800	102
7	Motala	Örebro	SUB	-6.600000	194.172879	92
8	Boden	Kiruna	SUB	-2.700000	433.841275	291
9	Sundsvall	Östersund	SUB	-2.300000	165.557147	166
10	Kalmar	Västervik	SUB	1.900000	108.463603	139
11	Halmstad	Helsingborg	SUB	-7.400000	174.313156	79
12	Umeå	Örnsköldsvik	SUB	13.299999	501.853856	111
13	Boden	Luleå	SUB	7.100000	50.642772	32
14	Lidköping	Skövde	SUB	6.100000	76.164332	49
15	Karlstad	Örebro	SUB	4.400000	104.013185	77
16	Halmstad	Varberg	SUB	11.399999	198.918661	65
17	Lidköping	Vänersborg	SUB	-2.700000	49.681729	60
18	Motala	Skövde	SUB	-7.100000	297.793737	118
19	Uddevalla	Vänersborg	SUB	1.200000	13.651827	21
20	Gävle	Sandviken	SUB	-3.500000	23.834965	25
21	Linköping	Motala	CORE	-110.900000	510.000000	51
22	Borås	Göteborg	SUB	-6.400000	131.078630	71
23	Falun	Västerås	SUB	4.600000	220.655492	128
24	Linköping	Västervik	SUB	-5.200000	167.094487	97
25	Sundsvall	Örnsköldsvik	SUB	-15.499999	711.515474	127
26	Jönköping	Motala	SUB	26.600000	954.098412	108
27	Vetlanda	Växjö	SUB	-10.600000	214.788040	72
28	Borås	Varberg	SUB	-14.600000	365.444485	84
29	Haparanda	Luleå	SUB	1.100000	58.031396	124
30	Luleå	Skellefteå	SUB	9.100000	451.386288	133
31	Stockholm	Uppsala	SUB	-29.400000	542.869621	69
32	Motala	Norrköping	CORE	-53.600000	790.000000	79
33	Kristianstad	Växjö	SUB	3.900000	172.119956	120
34	Motala	Vetlanda	SUB	-13.400000	674.669033	135
35	Gävle	Uppsala	SUB	26.400000	397.999170	60
36	Helsingborg	Malmö	SUB	-4.400000	74.666529	60
37	Borås	Jönköping	SUB	24.400000	582.941924	82
38	Malmö	Trelleborg	SUB	-3.600000	32.569445	34
39	Hudiksvall	Sundsvall	SUB	-18.799999	443.396458	81
40	Eskilstuna	Norrköping	SUB	9.400000	316.216164	102

core_cost = sum(edge_info[edge_info['Type'] == 'CORE']['Cost'])
subnet_cost = sum(edge_info[edge_info['Type'] == 'SUB']['Cost'])
total_cost = core_cost + subnet_cost

print(f'core   = {core_cost:>9.3f}')
print(f'subnet = {subnet_cost:>9.3f}')
print('------------------')
print(f'total  = {total_cost:>9.3f}')
assert abs(autonomax.model.getObjective().getValue() - total_cost) < 1e-6

core   =  1740.000
subnet = 11955.037
------------------
total  = 13695.037

visualization.show(pd.DataFrame(autonomax.city_info()), pd.DataFrame(autonomax.edge_info()))